#include "VertMagnDipole.h"
#include "mfAddFuncCalcTask.h"
#include "BesselFunctions.h"
#include "kernel/feObserver.h"
#include "kernel/feDefs.h"
#include "FastMath.h"
#include <math.h>
#include "Timer.h"

//---------------------------------------------------------------------------------------------------------------------------------------------
const InfoM  VertMagnDipole::InfosM[11] =
{
	{   1.000000e-007  ,  1.000e-003  ,  1.000e+1   ,  1.05  } ,
	{   1.000000e-006  ,  1.000e-004  ,  1.000e+1   ,  1.05  } ,
	{   1.000000e-005  ,  1.000e-005  ,  1.000e+1   ,  1.05  } ,
	{   1.000000e-004  ,  1.000e-006  ,  1.000e+1   ,  1.05  } ,
	{   1.000000e-003  ,  1.000e-06   ,  1.000e+1   ,  1.05  } ,
	{   1.000000e-002  ,  1.000e-05   ,  1.000      ,  1.05  } ,
	{   1.000000e-001  ,  7.000e-013  ,  1.000e-1   ,  1.05  } ,
	{   1.000000e+000  ,  7.000e-013  ,  1.000e-2   ,  1.05  } ,
	{   1.000000e+001  ,  5.000e-013  ,  1.000e-003 ,  1.05  } ,
	{   1.000000e+002  ,  1.000e-014  ,  1.000e-006 ,  1.05  } ,
	{   1.000000e+005  ,  1.000e-015  ,  1.000e-007 ,  1.05  }
};


//---------------------------------------------------------------------------------------------------------------------------------------------
VertMagnDipole::VertMagnDipole()
{
	MsCount = 0;
	MsArray = 0;

	Layers = 0;
	NumLayers = 0;

	DzDts = 0;

	TimesCount = 240;
	double tau = 1.1;
	double ftime = 1e-10;
	double step = ftime;
	double t = ftime;	

	//fill time array
	TimesArray = new double[TimesCount];
	for(int i=0; i<TimesCount; i++)
	{
		TimesArray[i] = t;
		step *= tau;
		t += step;
	}

	MaxAnalT = 0;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
void VertMagnDipole::TimeInterpolation(vector<double> &iFrom, double *oTo, double *iTimes, int iTimeCnt)
{
	//interpolation from main time mesh to extern mesh
	int it = 0;
	for(int i=1; i<TimesCount-1 && it < iTimeCnt;)
	{
		double t = iTimes[it];
		if(TimesArray[i] > t)
		{			
			double ts[3] = {TimesArray[i-1], TimesArray[i], TimesArray[i+1] };
			double fs[3] = { iFrom[i-1], iFrom[i], iFrom[i+1] };
			oTo[it] = feInterpolation::Quadratic(ts, fs, t);
			it++;
		}
		else
			i++;
	}
}


//---------------------------------------------------------------------------------------------------------------------------------------------
inline double VertMagnDipole::BesselInt(double iF1, double iF2, double iR, double iX0, double iX1)
{
	BesselFunctions& bf = BesselFunctions::GetInstance();


	double a = iF1/(iX0-iX1)-iF2/(iX0-iX1);
	double b = -iX1*iF1/(iX0-iX1)+iX0*iF2/(iX0-iX1);
	

	double j0x0 = bf.J0(iX0*iR);
	double j1x0 = bf.J1(iX0*iR);
	double s0x0 = bf.StruveH0(iX0*iR);
	double s1x0 = bf.StruveH1(iX0*iR);

	double j0x1 = bf.J0(iX1*iR);
	double j1x1 = bf.J1(iX1*iR);
	double s0x1 = bf.StruveH0(iX1*iR);
	double s1x1 = bf.StruveH1(iX1*iR);


	double res = -0.5*(2*a*iX0*j1x0 + 2*b*iR*iX0*j0x0 + b*PI*iR*iX0*s0x0*j1x0
		- b*PI*iR*iX0*s1x0*j0x0 - 2*a*iX1*j1x1 - 2*b*iR*iX1*j0x1 -
		b*PI*iR*iX1*s0x1*j1x1 + b*PI*iR*iX1*s1x1*j0x1)/iR;
	return res;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
inline double VertMagnDipole::IntFunction(int iMsIndex, int iTimeIndex)
{
	//done
	return MsArray[iMsIndex]*MsArray[iMsIndex]*DzDts[iTimeIndex][iMsIndex];
}


//---------------------------------------------------------------------------------------------------------------------------------------------
void VertMagnDipole::SetMedium(ConductiveLayer *iLayers, int iNumLayers)
{
	NumLayers = iNumLayers;
	Layers = iLayers;

	//get max time where analytic solvation is possible
	if(NumLayers == 1)
	{
		MaxAnalT = TimesArray[TimesCount-1]*2;
		return;
	}

	NeedRecalc = true;
// 	//iterating over infosM table
// 	for(int i=10; i>=0; i--)
// 	{
// 		//max height
// 		double h = 30/InfosM[i].M;
// 
// 		if(h < Layers[0].H)
// 		{
// 			if(InfosM[i].Tmin <= TimesArray[0])
// 				MaxAnalT = InfosM->Tmax;
// 		}
// 	}
}


//---------------------------------------------------------------------------------------------------------------------------------------------
bool VertMagnDipole::CreateDzDt(double iZ)
{
	if(!Layers)
		return false;


	//delete m
	if(MsArray)
		delete [] MsArray;
	MsCount = 0;

	//m-mesh params
	//7 per 5 degrees = 75 + 50x3 from 1 to 100
	MsCount = 100;

	double pointsperdegree = 11;
	double startdegree = -6;


	//allocate ms
	MsArray = new double[MsCount];
	
	//allocate dzdts
	if(!DzDts)
	{
		DzDts = new double *[TimesCount];
		for(int it=0; it<TimesCount; it++)
			DzDts[it] = new double[MsCount];
	}


	mfAddFuncCalcTask ztask;
	ZTaskHelper helper;
	
	double z = iZ < 0 ? 0 : iZ;
	
	//fill array
	double m = 1e-6;
	int iinf = 0;
	for(int im=0, i=0; im<MsCount; im++, i++)
	{		
		m = pow(10, startdegree+i/pointsperdegree);		

		//not const mesh params
		if(m >= 1 && pointsperdegree < 30) 
		{
			startdegree = 0; 
			i = 0;			
			pointsperdegree = 40;
			m = pow(10, startdegree+i/pointsperdegree);
		}


		MsArray[im] = m;

		//read task params from table
		while(InfosM[iinf].M < m) iinf++;
		double firstt = max(InfosM[iinf].Tmin, TimesArray[0]);
		double maxt = min(InfosM[iinf].Tmax, TimesArray[TimesCount-1]);
		double timecoef = InfosM[iinf].Tau;

		//get num layers over geometric progression sum
		int nl = (int)(log(1-maxt/firstt * (1-timecoef))/log(timecoef) + 0.5);
		
		helper.Init(Layers, m);

		//init task
 		if(!ztask.Init(Layers, NumLayers, timecoef, firstt, nl, m, &helper, true)) continue;
 
		//solve
		if(!ztask.Solve()) continue;


		//for z < 0
		double mult = iZ < 0 ? exp(m*iZ) : 1;


		//fill polyline for current m
		for(int it=0; it<TimesCount; it++)
			DzDts[it][im] = ztask.GetDuDt(z, TimesArray[it]) * mult;		
	}

	NeedRecalc = false;
	return true;
}


//--------------------------------------------------------------------------------------------------------------------------------------
bool VertMagnDipole::Br(double iR, double iZ, double *oResV, double *iTimes, int iTimeCount) 
{
	return false;
}


//--------------------------------------------------------------------------------------------------------------------------------------
bool VertMagnDipole::Bz(double iR, double iZ, double *oResV, double *iTimes, int iTimeCount) 
{
	return false;
}


//--------------------------------------------------------------------------------------------------------------------------------------
bool VertMagnDipole::Ephi(double iR, double iZ, double *oResV, double *iTimes, int iTimeCount) 
{
	return false;
}


//--------------------------------------------------------------------------------------------------------------------------------------
bool VertMagnDipole::DBrDt(double iR, double iZ, double *oResV, double *iTimes, int iTimeCount)
{
	return false;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
double VertMagnDipole::AnalDbzDt(double iT, double iR)
{
	const double sigma = Layers[0].Sigma;
	const double c1 = 9/(2*PI * sigma * pow(iR, 5));
	const double c2 = FastMath::Sqrt(2/PI);
	const double u = iR * FastMath::Sqrt(MU0 * sigma / (2*iT));
	
	double phi = 0;

	//simpson method
	const int n = 200;
	double h = u/n;


	phi = 1;
	for(int i=1; i<n; i++)
	{
		double x = i*h;
		double f = exp(-x*x/2);

		if(i%2 == 0)
			phi += 2 * f;
		else
			phi += 4 * f;
	}

	const double lastf = exp(-u*u/2);
	phi += lastf;
	phi *= c2 * h/3;


	return -c1 * (phi - c2 * u * (1 + u*u/3 + u*u*u*u/9) * lastf);
}



//--------------------------------------------------------------------------------------------------------------------------------------
bool VertMagnDipole::DBzDt(double iR, double iZ, double *oResV, double *iTimes, int iTimeCount)
{
	if(!Layers) return false;
	if(!DzDts || NeedRecalc) CreateDzDt(iZ);

	std::vector<double> dbzdt;
	dbzdt.reserve(TimesCount);


	//try to extract from cache
	if(BzDtCache.GetV(iR, 0.01*iR, dbzdt))
	{
		TimeInterpolation(dbzdt, oResV, iTimes, iTimeCount);
		return true;
	}


	//calc new
	double zero_t = Layers[0].H*Layers[0].H  * MU0 * Layers[0].Sigma;
	for(int it = 0; it<TimesCount; it++)
	{
		//analytics while in first layer
// 		if(TimesArray[it] < zero_t)
// 		{
// 			dbzdt.push_back(AnalDbzDt(TimesArray[it], iR));
// 			continue;
// 		}


		//numerical
		double s = 0;
		double f[2];
		double m[2];

		m[0] = MsArray[0];
		f[0] = IntFunction(0, it);
		

		for(int i=1; i<MsCount; i++)
		{
			//fill arrays
			f[1] = IntFunction(i, it);
			m[1] = MsArray[i];

			s += (m[1]-m[0])*0.5*(f[0]*BesselFunctions::GetInstance().J0(iR*m[0]) + f[1]*BesselFunctions::GetInstance().J0(iR*m[1]));//BesselInt(f[0], f[1], iR, m[0], m[1]);	

			//update f[0]
			f[0] = f[1];
			m[0] = m[1];			
		}


		//push solvation
		s *= -MU0/(4*PI);
		dbzdt.push_back(s);
	}

	BzDtCache.PushV(iR, dbzdt);
	TimeInterpolation(dbzdt, oResV, iTimes, iTimeCount);

	return true;
}


//--------------------------------------------------------------------------------------------------------------------------------------
bool VertMagnDipole::DEphiDt(double iR, double iZ, double *oResV, double *iTimes, int iTimeCount) 
{
	return false;
}


//--------------------------------------------------------------------------------------------------------------------------------------
VertMagnDipole::~VertMagnDipole()
{
	if(TimesArray) delete [] TimesArray;
}